function F = solve_eqm_inner(param, const, integ)
%% Initial Equilibrium

    % normalization
    P_s     = param.P_s;
    e       = param.e;     
    X       = param.X;

    % read in parameters
%     kappa   = param.kappa;
    alpha_m = param.alpha_m; 
    alpha_s = param.alpha_s;
%     beta_v  = param.beta_v; 
%     beta_m  = param.beta_m;    
    sigma   = param.sigma; 
    sigma_s = param.sigma_s;
%     f_vH    = param.f_vH;
%     f_vF    = param.f_vF;
    mu_E    = param.mu_E;
    sigma_E = param.sigma_E;   
    w       = param.w;
    Q       = param.Q;
    Nf      = param.Nf;
    
    % read in constants
    gamma   = const.gamma;
    k0      = const.k0;
%     C_m     = const.C_m;
%     C_v0    = const.C_v0;
%     C_x0    = const.C_x0;   
%     phi_v   = const.phi_v;
    phi_y   = const.phi_y;
    phi_s   = const.phi_s;
    D_F     = const.D_F;
    lnzQ    = const.lnzQ;
        
    % read in all the integral from previous round
%     E_zm0 = integ.E_zm0;
%     E_zm1 = integ.E_zm1;
%     E_zv0 = integ.E_zv0;
%     E_zv1 = integ.E_zv1;
    E_zx0 = integ.E_zx0;
    E_zx1 = integ.E_zx1;
        
    % adjust for level in case integrals too large
    rescale = max(1,10^floor(log10(max(E_zx1))-1));
    E_zx0 = E_zx0/rescale;
    E_zx1 = E_zx1/rescale;
%     E_zv0 = E_zv0/(rescale^(1/beta_v));
%     E_zv1 = E_zv1/(rescale^(1/beta_v));
%     E_zm0 = E_zm0/(rescale^(1/beta_m));        
%     E_zm1 = E_zm1/(rescale^(1/beta_m));      

    
    newmbar = 1.072121591267341e-04;


%% Solve the equilibrium    

    % Setting
    tol = 1e-12;
    maxIter = 1000;

    % Initial guess
    D_Hnew = ones(1,Q); 
    c_new = ones(1,Q);

    % Iteration
    err = 1;
    iter = 0;
    while err >= tol && iter < maxIter

        % new guess
        D_H = D_Hnew;
        c = c_new;

        % D(q,E)
        D0 = D_H;
        D1 = D_H + e^sigma.*D_F;        
        
        % Pi(q)
        C_x = (sigma/(sigma-1)*w.^(1-alpha_m-alpha_s)).^(1-sigma);
        Pi0 = D0.*(c.^alpha_m*P_s^alpha_s).^(1-sigma).*C_x;
        Pi1 = D1.*(c.^alpha_m*P_s^alpha_s).^(1-sigma).*C_x;
 
        % X(q,E)
        X0 = Pi0.*E_zx0;
        X1 = Pi1.*E_zx1;
        X_m = alpha_m*(sigma-1)/sigma*(X0+X1);

        % P(q)^(1-sigma)
        P0sig = X0./D0;
        P1sig = X1./D1;
        Psig = X0./D0 + X1./D1;

        % D_m(q)
        TempD = c.^(sigma-1).*X_m;
        TempD(isnan(TempD))=0;
        D_m = newmbar*sum(phi_y.*repmat((TempD)',1,Q),1);

        % spending on domestic services
        B1 = sum(e^sigma*D_F./D1.*X1); %export       
        X_sH = (1-(sigma-1)*alpha_m/sigma)*X-min(B1,0.7);       
        
        % D_s(q)
        D_s = phi_s.*X_sH/sum(phi_s.*Psig);

        % D(q) (new guess)
        D_Hnew = D_m+D_s;

        % c(q) (new guess)
        c_new = (newmbar.*sum(phi_y.*repmat(Psig,Q,1),2)').^(1/(1-sigma));   
        
        % update
        D_Hnew = D_H + 0.2*(D_Hnew - D_H);
        c_new = c + 0.2*(c_new - c);
        iter = iter+1;      

        % check convergence
        err1 = max(abs(D_Hnew-D_H)./D_H);
        err2 = max(abs(c_new-c)./c);
        err = err1 + err2;

    end
    
    % Warning
    if iter == maxIter
        msg = 'User Warning: Inner Loop Maximum Iterations Reached.';
        warning(msg)
    end 

    % Adjust back
    Pi0 = Pi0/rescale;
    Pi1 = Pi1/rescale; 
    
    % Save the Export Decision
    ExportCutoffQ = (k0*lnzQ-log(sigma*gamma)+log(kron(ones(Nf,1),Pi1)-kron(ones(Nf,1),Pi0))-mu_E)/sigma_E;
    ExportProbQ = normcdf(ExportCutoffQ);        


    
%% Return

    eqm.Pi0 = Pi0;
    eqm.Pi1 = Pi1;
    eqm.D_H = D_H;
    eqm.D_F = D_F;
    eqm.D0 = D0;
    eqm.D1 = D1;
    eqm.D_s = D_s;
    eqm.D_m = D_m;
    eqm.c = c;
    eqm.Psig = Psig;
    eqm.P1sig = P1sig;
    eqm.P0sig = P0sig;
    eqm.P = Psig.^(1/(1-sigma));
    eqm.X_m = X_m;
    eqm.X0 = X0;
    eqm.X1 = X1;
%     eqm.V = V;
%     eqm.Vbar = Vbar;
%     eqm.V0 = V0;
%     eqm.V1 = V1;
%     eqm.M = M;
%     eqm.M0 = M0;
%     eqm.M1 = M1;
%     eqm.xi = xi;
%     eqm.theta_v = theta_v;
%     eqm.theta_m = theta_m;  
%     eqm.rv_ads = rv_ads;
%     eqm.rv_rev = rv_rev;
%     eqm.C_v1 = C_v1;    
    eqm.B1 = B1;
    eqm.P_s = P_s;
    eqm.e = e;
    eqm.X_sH = X_sH;
    eqm.X = sum(X0+X1);
%     eqm.P_sF = P_sF;
%     eqm.P_sH = P_sH;
%     eqm.mbar = mbar;  
%     eqm.f = f;
    eqm.ExportCutoffQ = ExportCutoffQ;
    eqm.ExportProbQ = ExportProbQ;
    eqm.iter_inner = iter;
    eqm.err_inner = err;
    
    eqm.newmbar = newmbar;
    
    F = eqm;

    
end